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^ ; ABSTRACT 

. Member galaxies within galaxy clusters nowadays can be routinely identified in cosmologi- 

Qh' cal, hydrodynamical simulations using methods based on identifying self bound, locally over 

q \ dense substructures. However, distinguishing the central galaxy from the stellar diffuse corn- 

el ■ ponent within clusters is notoriously difficult, and in the center it is not even clear if two 

' distinct stellar populations exist. Here, after subtracting all member galaxies, we use the ve- 

& . locity distribution of the remaining stars and detect two dynamically, well-distinct stellar com- 

ponents within simulated galaxy clusters. These differences in the dynamics can be used to 
apply an un-binding procedure which leads to a spatial separation of the two components 
into a cD and a diffuse stellar component (DSC). Applying our new algorithm to a cosmo- 

■ logical, hydrodynamical simulation we find that - in line with previous studies - these two 
| components have clearly distinguished spatial and velocity distributions as well as different 

star formation histories. We show that the DSC fraction - which can broadly be associated 
with the observed intra cluster light - does not depend on the virial mass of the galaxy clus- 
ter and is much more sensitive to the formation history of the cluster. We conclude that the 
separation of the cD and the DSC in simulations, based on our dynamical criteria, is more 
' physically motivated than current methods which depend on implicit assumptions on a length 

scale associated with the cD galaxy and therefore represent a step forward in understanding 
the different stellar components within galaxy clusters. Our results also show the importance 
of analyzing the dynamics of the DSC to characterize its properties and understand its origin. 
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1 INTRODUCTION 

The existence of a diffuse stellar component in groups and 
clusters of galaxies is now well established. A diffuse intra- 
cluster light (ICL) has been observed both in local universe 
dFeldmeier eTal] |2004 iMihos et all 12001 lArnaboldi etail |2004| ; 
iGerhard et al. 2005; Dohertv et a l. 2009) and at inte r mediate red- 



shift i Gonzalez et alj l2000l: iFeldmeier etail |2004 Izibetti et al] 
120051 : iKrick & Bernsteinll2007l)~ Such ICL component is centrally 
conce ntrated, and typically amounts bet ween ss 10 jZibetti et al.1 
l2005h and « 35 % ^Gonzalez et al.l2007t) of the total stellar mass in 
clusters. There are some indications pointing towards a dependence 
of the ICL fraction from the mas s of the clusters, going from ~ 2 
% at the scale of loose groups (Castro-Rodriguez et al. 2003h to 
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« 5-10 % of Virgo cluster iArnaboldi et al.l20ollFeldmeier et al] 
2004lMihos et alfcOOlrCastro-Rodriguez et alj|2009h up to 10 - 
20 % or higher in the most massive clusters ^Gonzalez et al.|200C ; 



Feldmeier et al. 



20021 ; iGal-Yam et aI1l2003l : Feldmeie r etai]|2004 



Krick et alteOOaT but other studies found instead no significant de 



pendence on the cluster richness jZibetti et alj2005h 

ICL is usually detected using surface photometry in various 
bands, s ometimes stacking together observations of many clusters 
(see e.g. Izibetti et al.| [2005). In this case, light coming from galax- 
ies is mas ked away, often using a fixed limit in surface brightness 
(e.g., as in Feldmeier et al. 2004). Single ICL stars can also be de- 
tected (Dur rell et alj|2002l) , expecially by observ i ng Intra-Cluster 
Plane tary Nebulae (ICPN) ( lArnaboldi et alJI 19961 : IFeldmeier et al] 
12004b. which can be observed up to distances as large as 100 Mpc 
dGerhard et"al] 12005). On the one hand, the observational effort 
needed to detect single stars makes it difficult to use them to as- 
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sess general properties of the ICL distribution, like the ICL frac- 
tion in clusters. On the other hand, the observation of single stars 
make it possible to gather information on the kinematic of such a 
component. In particular, extended halos of bright ellipticals often 
overlap spatially with stars which are free-floating in the cluster 
gravitational potential. Effo rts have been made to disentangle these 
two components (see e.g. iDohertv et al] |2009), but they are often 
referred together as the ICL. In the following, we will use ICL to 
indicate the observed diffuse stellar population, whose dynamical 
properties are difficult to obtain, and DSC for the corresponding 
stellar population produced in numerical simulation, whose dynam- 
ics is known. 

From a theoretical point of view, a diffuse stellar component 
(DSC) has been st udied in a cosmolo gical Dark-Matter-only nu- 
merical simulation (jNapolitano et alT20 03). Later, for the first time 
iMurante et al.l d2004h (M04 hereafter) studied the properties of the 
DSC in a cosmological simulation which included hydrodynam- 
ics and various astrophysical processes s uch as radiative coolin g of 
gas and (self-consistent) star formation. IWillman et alj J20041) and 
ISommer-Larsen et al] d2005h analyzed a DSC in their resimulation 
of single galaxy clusters. 

M04 analyzed a cosmological simulation of a box of 192 h~ x 
Mpc, containing w 100 galaxy cluster. Their main results were 
that: (i) a DSC is clearly seen when a Sersic fit is performed on 
the bound, unbound and total two-dimensional radial density pro- 
files of the stellar population in clusters; (ii) the DSC component is 
more centrally concentrated than the stars bound in galaxies; (iii) 
on average, stars in the DSC are older than th ose in galaxies; ( iv) the 
DSC fraction increase with the cluster mass. IWillman et all d2004h 
used one cluster re-simulation, having a mass of 1.2 x 10 15 M Q , 
and confirmed results (i) of M04. They found a DSC fraction com- 
patible with M04 clus t ers in the same mass range, if slightly lower. 
ISommer-Larsen et alJ ( Eoolh used one Virgo-like and one Coma- 
like cluster re-simulation and confirmed results (i), (ii) and (iii) of 
M04. They found similar DSC fraction, with their Coma-like clus- 
ter having an higher fraction of DSC stars. Their average stellar age 
in DSC was older than that found by M04. 

Several mechanisms have been investigated to explain the 
formation and evolution of a DSC: stripping and disruption of 
gala xies as they pass through the central regions of relaxed clus- 
ters l lBvrd & Valtonenlll990l;lGnedinir2 003 ) ; stripp ing of stars from 
galaxies during the initial formation of clusters dMerrittl 1984b : cre- 
ation of stellar haloes in galaxy group s, that later fall into mas- 
sive c lusters and then become unbound (Mihos 2004; Rudick et al.l 
2006); and stripping of stars during high-speed g alaxy encoun- 
ters in the cluster environment dMoore et alj|l996t) . IMurante et al.l 
d2007l) (M07 hereafter) showed that the DSC is produced during 
mergers in the formation history of the BCGs and of other massive 
galaxies, and that it grows st eadily since redshift z = 1, with no 
preferred epoch of formation. iMonaco et al.l d2006l) . using a semi- 
analytical model of galaxy formation in clusters, showed that, in 
the hierarchical model of structure formation, the BCG merging 
activity alone between z = 1 and z = is so intense that it is 
not possible to simultaneously fit the bright end of the galaxy lumi- 
nosity function at both redshifts without assuming that a significant 
fraction of stars goes into the DSC during such mergers. We point 
out that observational indications have be en found that such p rocess 
is in fact in action in the Coma clu ster dGerhard et alj|20071) . Nu- 
merically, IStanghellini et al.l d2006l) studied isolated, star-only dry 
merger of elliptical galaxies and showed how in that case, up to 21 
the initial total stellar mass can become unbind in the process. 

One major problem with the analysis of DSC in numerical 



simulations is its ve ry identification. One possibility, used e.g. by 
Rudick et al] J2006h . is to build simulated surface brightness maps 
and detect an ICL component with a procedure similar to that used 
with observational data. This has the advantage of being somewhat 
easier to compare with observation, at the cost of losing the intrin- 
sic advantage of simulations, that is, the knowledge of the dynamics 
of each star particle. The other possibility is to attempt a dynami- 
cal distinction between stars bound to galaxies and stars which are 
free-floating in the g rou p or cluster potential. This i s what M04, 
Will man et ail d2004b and lSommer -Larsen ct ai] d2005h did; the dis- 
advantage is that such a distinction may or may not coincide with 
the observational definition of ICL. However, for doing this, it is 
necessary to identify structures in the simulation (groups or clus- 
ters) and substructures within them (the galaxies). Up to date, no 
unambiguous way to perform such identification has been formu- 
lated. 

The identification of structures is a relatively easy task. The 
most commonly used algorithms ar e based upon two methods: the 



Frien ds-of-Friends (FoF) approach dDavisetalJll984 



, Frenk et al 

Il988h and the Spherical Overdensity one (SO) dLacev & Cola 
1 19941) . In the former, all particles closer than a given distance I 
(usually expressed in units of the mean inter-particle distance, typ- 
ically I w 0.2 is chosen) are grouped together; then, all particles 
closer than I to such particles are included in the group, and so on. 
In the latter, one identifies density peaks in the particle distribution, 
and puts spheres around them, choosing their radius so that they 
enclose a given overdensity. 

The difficulty arises when ™£>structures inside DM halos must 
be found. The density contrast of substructures against their parent 
halo density is lower than the density contrast of structures against 
the background, and this makes their identification more problem- 
atic. Even worse, the gravitational potential of a structure has a 
clear zero-point, given by the background density of the universe. 
So, it is well defined which particles are bound to the structure and 
what are unbound. Inside a structure, instead, such a definition is 
not clear: the background density of the structure varies inside them 
with the distance from the center, the dynamical state of the cluster, 
and the typical extend of the substructure itself. Therefore, telling 
that a particle is gravitationally bound to a given substructures in- 
volves a certain degree of arbitrariness, and an operational defini- 
tion must be used. 

A number of different subhalo finders now exist in the 
literature. Some are based on pure geometrical measures, oth- 
ers construct a candidate structure and apply an unbinding cri- 
teria, e.g. comparing the kinetic energy of a particle to its 
potent ial energy. A n on-exhaustive list includes: hierarchical 
FOF dKlypin et alj Il999l). Bound Density Maxima dklypin et al 



1999b . DENMAX (Bertschin ger & Gelb 1991; Gel b & Bertsch inger 
1994b . SKID dWeinberg et al.l 1 19971: IStadell 1200 lb. local den- 
sity percolation (vanKampen 1995), HOP (Eisenstein& Hull 
1 1998b . SUBFI ND dspringeletallBjOlh . MHF IGill et alj l20O4b . 



ADA PTAHOP dAubert et alj|2004l) . V OBOZ dNevrinck et al 



2005b . 



PSB dKim & Parkl 120061) 7 6D FOF iDiemand et alj d2 006), HSF 
(Maci eiewski et al.l2009b . AHF dKnollmann & Knebel2 009). Some 
of these algorithms are based on the FoF scheme, and vary the 
I parameter inside structures or apply the scheme to the full 6D 
phase space (hierarchical FOF, local density percolation, 6D FOF); 
some use the SO scheme, and apply an unbinding criterion (BDM, 
MHF, AHF); some use FoF to link particles to density maxima, 
and apply an unbinding (DENMAX, SKID); some introduce geo- 
metrical constraints, together with various unbinding procedure and 
phase space analysis (HOP, SUBFIND, ADAPTAHOP, VOBOS, PSB, 
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HFS). Comparisons between the performance of the different algo- 
rithms in detecting substructures usually show a reasonable agree- 
ment between the different m et hods as can be found for e xample in 
Knollmann & Knebej J2009I) or lMacieiewski et alj J2009h . Already 
Spring el et al.l d200lh compared SUBFIND and FOF, showing that 



the former performs better in those cases where nearby two halos 
are connected by a filaments. FOF usually joins together the two 
objects, while SUBFIND can disentangle them. Such a behaviour is 
however common to most of the subhalo finder listed above. 

As far as the DSC detection is concerned, a further question 
arises. The DSC is composed by stars bound to the gravitational 
potential of the cluster and not to any galaxy, and it is centrally 
concentrated. Many of its general properties, therefore, are primar- 
ily determined by the stars at the center of the cluster. It is therefore 
important to correctly distinguish DSC star particles from cD star 
particles. But such a distinction is often somewhat arbitrary, as in 
SKID, where the unbinding procedure is based upon the value of 
a parameter that must be given externally and whose effect must 
be tested a-posteriori (see M07). Other substructure finders, like 
SUBFIND, identifies all substructures but the main one and link all 
remaining star, gas and DM particles to the main subhalo, thus mix- 
ing cD and DSC stars. It is not even clear if a distinction between 
cD and its extended stellar halo is physically motivated or if it is all 
arbitrary. 

In this work, we show that such two components can be dis- 
tinguished in the velocity space. We use SUBFIND to identify sub- 
structures and disentangle galaxies from the cD + DSC stars. We 
show that the velocity distribution of the latters cannot be fit with 
a single Maxwellian distribution, while it is well fitted by a double 
Maxwellian (the sum of two Maxwellian distribution with different 
velocity dispersions). 

We thus modify SUBFIND to perform an unbinding procedure 
on the cD + DSC stars. We use the gravitational potential given 
by the matter contained in a sphere, centered on the halo center, 
and we find a radius which divides the star particles in two pop- 
ulations whose velocity distribution is separately fit by a single 
Maxwellian. We vary the radius until the two resulting velocity 
dispersions correspond to the two velocity dispersion given by the 
double Maxwellian fit of the velocity distribution of the whole cD 
+ DSC star population. We show that this procedure is able to dis- 
entangle the central galaxy from its extended stellar halo. We thus 
characterize the DSC through its phase-space distribution. 

We then use our modified SUBFIND algorithm to repeat the 
M04 analysis of the DSC applied to a constrained cosmological 
simulation of the local universe. We also use two re-simulations of 
galaxy clusters, to test the behavior of our modified SUBFIND when 
we i ncrease resolutio n. 

IPuchwein et alj |2010) used different ways to define the DSC 
in their simulations, including the one we present here. They found 
a broad qualitative agreement among the various methods. Here, 
we will focus on the comparison of our modified SUBFIND algo- 
rithm with SKID on our constrained cosmological simulation of the 
local universe, and show that they are in broad agreement, but dy- 
namically separating the cD and the DSC, a dependance of the DSC 
fraction on cluster mass is not found. 

The plan of the paper is as follows. Section [2] describes our 
simulations. Section|3]shows the existence of two dynamically dis- 
tinct stellar components, even at the center of galaxy clusters. Sec- 
tion [4] presents our scheme for disentangling such components. In 
Section [5] we apply our scheme to our local universe simulation, 
also comparing the results we obtain on the DSC properties with 
those obtained using SKID. In Section [6] we give our conclusions. 



In the Appendix we show a detailed comparison of the performance 
of SKID and SUBFIND for two galaxy clusters. 



2 GENERAL SETUP 
2.1 Simulations 

The results presented in this paper have been obtained using the 
final output of a cosmological hydrodynamical simulation of the 
lo cal universe. Our in itial conditions are similar to those adopted 
bv lMathis et alj (12002') in their study (based on a pure N-body sim- 
ulation) of structure formation in the local universe. The galaxy 
distribution in the IRAS 1.2-Jy galaxy survey is first Gaussian 
smoothed on a scale of 7 Mpc and then linearly evol ved back in 
time u p to z — 50 following the method proposed by Kolatt et al. 
( 1996) . The resulting fie ld is then used as a Gaussian constraint 

1 Hoffman & Ribaklll99lh for an otherwise random realization of a 
flat ACDM model, for which we assume a present matter density 
parameter flo-m = 0.3, a Hubble constant Ho = 70 km/s/Mpc 
and a r.m.s. density fluctuation eg = 0.9. The volume that is con- 
strained by the observational data covers a sphere of radius ~ 80 
Mpc/h, centered on the Milky Way. This region is sampled with 
more than 50 million high-resolution dark matter particles and is 
embedded in a periodic box of ~ 240 Mpc/h on a side. The re- 
gion outside the constrained volume is filled with nearly 7 million 
low-resolution dark matter particles, allowing a good coverage of 
long-range gravitational tidal forces. 

The statistical analysis made by Ma this et alj j2002) showed 
that, using this technique, the simulated z — matter distribution 
provides a good match the large-scale structure observed in the lo- 
cal universe. Moreover, many of the most prominent nearby galaxy 
clusters like Virgo, Coma, Pisces-Perseus and Hydra-Centaurus, 
can be identified directly with halos in the simulation. Therefore 
the galaxy clusters in this simulations are an ideal target to study 
the diffuse stellar component, since observations are mainly bound 
to study objects in the local universe. 

Unlike in the original simulation made bv lMathis et aljj2002l) . 
where only the dark matter component is present, here we want to 
follow also the gas and stellar component. For this reason we ex- 
tended the initial conditions by splitting the original high-resolution 
dark matter particles into gas and dark matter particles having 
masses of 0.48 x 10 9 M Q /ftand3.1 x 10 9 M©//i, respectively; 
this corresponds to a cosmological baryon fraction of 13 per cent. 
The total number of particles within the simulation is then slightly 
more than 108 million and the most massive clusters is resolved by 
almost one million particles. The gravitational force resolution (i.e. 
the comoving softening length) of the simulation has been fixed to 
be 7 kpc/h (Plummer-equivalent), fixed in physical units from z=0 
to z=5 and then kept fixed in comoving units at higher redshift. 

To get a handle on possible effects of resolution on our mod- 
ified SUBFIND identification scheme for the cD and the DSC, we 
also analyzed two re-simulations of a massive galaxy clusters at 
higher resolution with exactly the same physics included. In the 
high resolution run (labeled as 3x in the following) the individ- 
ual particle masses are a factor of 3 times smaller than in our 
local universe simulation, and we scaled the softening accord- 
ingly to 5 kpc/h. Besides the most massive mail halo (M„; r m 

2 x 10 Mq/K), this re-simulation hosts several less massive ha- 
los (M V i r « 1 x 10 14 Mq /h), which we analyzed as well. In the 
extreme high resolution run (labeled as 18x in the following) the 
individual particle masses are 18 times smaller and the softening is 
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Figure 1. Visualization of the stellar component using the the ray-tracing software SPLOTCH. The color composition reflect the three dimensional velocity 
field (see text for details), whereas the intensity reflects the stellar density. The upper panel shows all stars in the high resolution sphere. The lower left panel 
shows all stars in galaxies (including the cD) of the most massive cluster, whereas the lower right panel shows all stars not bound to individual galaxies. 
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accordingly set to 2.5 kpc/h. Here, due to additional optimizations 
to the high resolution region of the re-simulation, only the main 
halo is suited to be analyzed. This cluster is resolved by more than 
respectively 4 and 25 million particles within the virial radius in the 
3x and the 18.x simulation. 

Our simul ations have been carried out with GADGET-2 
dSpringelllijoll). The code uses an en tropy-conserving formula- 
tion of SPH dSpringel & He rnquist 2002), and allows a treatment 
of radiative cooling, heating by a UV background, and star for- 
mation and feedback processes. The latter is based on a sub- 
resolutio n model for the multiphase structure of the interstellar 
medium (Springel & Hernquist 2003). The code can also follow the 
pattern of metal production from the past history of cosmic star for- 
mation dTornatore et alj2004ll2007h . This is done by computing the 
contributions from both Type-II and Type-la supernovae and energy 
feedback and metals are released gradually in time, accordingly to 
the appropriate lifetimes of the different stellar populations. This 
treatment also includes, in a self-consistent way, the dependence 
of the gas cooling on the local metalic ity. The feedback scheme 
assumes a Salpeter IMF dSalpeteilll955ri . and its parameters have 
been fixed to get a wind velocity of « 480 km/s. 

We note that our present local universe simulation has ad- 
ditional physics with respect to M04 work. Their mass resolu- 
tion was 4.62 x 10 9 /i _1 M for a Dark Matter (DM) particle 
and 6.93 x 10 8 /i _1 Mq for a gas particles. The force resolu- 
tion was 7.5k" 1 kpc (Plummer-equivalent softening length). While 
force and mass resolution are similar to those of our local universe 
simulation, M04 didn't include chemical enrichment and metal- 
dependent cooling. 



2.2 Substructure detection 

Substructures within halos are usually defined as locally over- 
dense, self-bound particle groups identified within a larger parent 
halo. In our analysis, the identification of these substructures is 
perfor med by applying a modified version of the SUBFIND algo- 
rithm ( Dolag et al. 2009), which, in contrast to its original version 
dSpringel et alj 2001 ) works with different particle species (e.g. 
dark matter, gas and star particle s ). Details on the a lgorithm can 
be found in dSpringel et al]|200ll: fPolag et alj|2009h . In brief, as 
a first step, we employ a standard friends-of-friends (FoF) algo- 
rithm to identify the parent halos, with a linking lenght = 0.17 
times the mean interparticle distance. Within each FoF group we 
then estimate the density of each particle species by adaptive 
kernel-interpolation, using the standard SPH approach with a given 
number of neighboring particles. Using an excursion set approach 
where a global density threshold is progressively lowered, we find 
locally over-dense regions within the resulting density field, which 
form a set of substructure candidates. The outer 'edge' of the sub- 
structure candidate is determined by a density contour that passes 
through a saddle point of the density field; here the substructure 
candidate joins onto the background structure. In a final step, all 
substructure candidates are subjected to a gravitational unbinding 
procedure where only the self-bound part is retained. We use all the 
particles initially belonging to the substructure candidate for eval- 
uating the gravitational potential. For the gas particles we take also 
the internal thermal energy into account in the gravitational unbind- 
ing procedure. If the number of bound particles left is larger than a 
20 particles, we register the substructure as genuine sub-halo. We 
define the stellar component of such sub-halo as a cluster member 
galaxy. 

Although this works well for identification of satellite galax- 



ies, all star particles not bound to any satellite galaxy are associated 
with the main halo stellar component, as the particles are always 
bound to the cluster potential. Therefore, this procedure does not 
split such a stellar population into the stars belonging to the central 
galaxy (cD) and to the stellar diffuse componen t (DSC). 

Previous work (M04, M07) based on SKID dStadefeOOlT) intro- 
duced an empirical distance scale to limit the part of the main halo 
used to calculate the potential and then distinguish between the cD 
and DSC by applying an unbinding criteria to this truncated poten- 
tial. In section|3]we will show that the main halo stellar component 
has two components with two distinct dynamics. In section [4] we 
show that this can be used to modify the unbinding procedure to 
replace empirical distance scale by an inferred one, which optimize 
the splitting of the two components according to their dynamics. 

To give an impression of the dynamics of the stellar com- 
ponent in our sim ulations we modified the SPLOTCH package 
(Dolag et al. 2008), a Ray-tracing visualization tool for SPH sim- 
ulations. Instead of mapping a scalar value to a color table, we 
mapped the three velocity components directly into the RGB col- 
ors. To obtain a logarithmic like scaling but preserve the sign 
of the individual velocity components we used a mapping based 
on the inverse of the hyperbolic sin function, e.g. < r,g,b >oc 
asinh(u x ,y,z)- We calculated a smoothed density field based on all 
star particles using a SPH like kernel function. The resulting Ray- 
tracing images, are shown in Figure [T] The upper panel shows all 
stars in the high resolution sphere. Same color indicates the same 
3 dimensional velocity direction, therefore the overall color pat- 
tern nicely reveals the large scale velocity structure within the cos- 
mological simulation. The larger clusters however appear mostly 
white meaning that the internal velocity dispersion overcomes the 
ordered, large scale velocity pattern and therefore the color satu- 
rates into white. The lower left panel shows all stars in galaxies 
(including the cD) of the most massive cluster (the one with the 
yellowish halo also seen close to the center in the upper panel). 
This cluster has 238 satellite galaxies above our detection thresh- 
old. The different colors of the galaxies reflect their different orbital 
phase, whereas the yellow color of the cD and some other satellite 
galaxies reflects the large scale velocity field as seen in the upper 
panel. As expected, besides the mean velocity of the cluster, there 
is no large scale pattern visible in the galaxy orbits. Also the indi- 
vidual galaxies show a uniform color indicating that the individual 
velocity dispersion of the stars within the galaxies is much smaller 
than their orbital motion within the galaxy cluster. The lower right 
panel shows all stars not bound to individual galaxies. The small 
scale structures remaining in the DSC are either galaxies below our 
detection threshold, or the outer envelopes of large cluster mem- 
bers, which are not bound to the galaxies. Notice the pink, ring like 
structure to the upper left, which can be associated to the group of 
infalling galaxies visible at the same position with the same color in 
the left panel. However, most of the DSC appears homogeneously 
in white and yellow colors, indicating a mean velocity similar to 
the cD and to the large scale velocity field, but a larger velocity 
dispersion. 



3 DISTINGUISHING DIFFERENT DYNAMICAL 
STELLAR COMPONENTS 

3.1 The different stellar components 

As shown in M07, the DSC within 0.5 R V i r forms mainly by merg- 
ing of satellite galaxies with the cD galaxy. Therefore one could 
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its bi-modality as it can be well fitted by the superposition of two 
Maxwellian distributions 



[km/s] 



Figure 2. Velocity histogram (black histogram) of the main halo stellar 
component and a double Maxwellian fit to it (black line). The red and the 
blue line are showing the individual Maxwellian distributions of the two 
components. The gray, dashed line is our best fit with a single Maxwellian. 




M <Ir [Msol/h] 



Figure 3. Velocity dispersion of the three stellar components (DSC, cD and 
galaxies) as a function of the virial mass of the galaxy cluster for the sample 
of 44 clusters with at least 20 satellite galaxies. The velocity dispersion of 
the DSC (blue symbols) is only slightly smaller than the one of the galaxies 
(black symbols), whereas the velocity dispersion of the cD stars (red sym- 
bols) is significantly smaller. The dashed lines show the one sigma scatter of 
our best fit relations. Velocity dispersion shows here are obtained using the 
double Maxwellian fitfT] The results from the high resolution simulations 
3x and I8x are shown with triangles and stars, respectively. 



expect that the DSC still reflects the dynamics of the satellite galax- 
ies and thus might have a velocity dispersion close to that of the 
satellite galaxies. On the other hand, the stars belonging to the cD 
galaxy will have a much smaller velocity dispersion due to the re- 
laxation and merging processes, which should get rid of the excess 
orbital energy in order to form the galaxy. Indeed, the main halo 
stellar component shows clear signs of a bimodal distribution. Fig- 
uref2]shows the histogram of the modulus of velocities ("velocities" 
hereafter) of the stellar component of the main halo for one of the 
more massive cluster of the simulation. This distribution reveals 



N(v) = kiv^exp ( ) + fc 2 t>fexp ( --§■ 
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For comparison, we also show in Figure [2] the single Maxwellian 
distribution which best fits our data. It is clear that while a single 
Maxwellian is a poor fit to the velocity distribution of the stellar 
component of our main halo, a double Maxwellian is a very good 
one. 

In itself, this does not allow to identify which star particle be- 
longs to the cD and which to the DSC but it already allows to draw 
conclusions based on the statistical properties of the two distribu- 
tions. Therefore we selected from our simulation all clusters having 
at least 20 identified satellite galaxies, ending up with a sample of 
44 galaxy clusters. For each of them we fitted a double Maxwellian 
to the velocity distribution of the main halo stellar component. Note 
that in all cases the fits are quite good, and introducing a third 
Maxwellian component does not improve them, as the amount of 
star particle associated to such third component always stays at the 
percent level, and just reflects structures like the outer envelope of 
large, in-falling groups. 

We hypothesize that the two Maxwellian distribution of our 
fits correspond to two distinct stellar components, namely the cD, 
associated with the smaller velocity dispersion, and the DSC hav- 
ing a larger one. Figure [3] shows the velocity dispersion of such 
two stellar components, together with the velocity dispersion of 
the galaxy population of each cluster, as a function of the virial 
mass of the galaxy cluster for our sample. The velocity dispersion 
of the DSC is only slightly smaller than the one of the galaxies: 
another hint on the origin of the DSC component from merging 
satellite galaxies. The slightly smaller velocity dispersion could 
indicate that either some momentum transfer occurred during the 
merger event with the cD, or the galaxies which contributed to the 
build up of the DSC till redshift zero have a slightly lower veloc- 
ity dispersion than the average galaxy population. The latter could 
be the case, since such galaxies might be mostly on radial orbits 
which lead to an earlier destruction by interactions with the cD than 
the (still) remaining galaxies within the cluster. However, since the 
DSC is more centrally concentrated than galaxies, they could sam- 
ple different region of clusters and thus have slightly different dy- 
namics. Figure [4] shows the radial velocity dispersion profile stack- 
ing all 44 clusters of our sample. From this analysis it appears that 
the dynamical differences are evidently the same at all radii. This 
demonstrates that such differences are not due to the fact that dif- 
ferent components sample different cluster regions. 

If the creation of the DSC would be a continuous process, one 
would expect a small scatter in the relation between virial mass of 
the cluster and velocity dispersions for the DSC component, spe- 
cially compared to that for galaxies, as they are much less numer- 
ous. However, from Figure [3] this does not seem to be the case, 
where we find a scatter of ~9% for the galaxies, «9% for DSC and 
~18% for the cD, again indicating that the production of the DSC 
originates in individual, violent events such as mergers with the cD. 
The velocity dispersion of the cD component is significantly lower, 
and shows an even larger scatter. In Figure [3] solid lines represent 
the best fit to 
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Figure 4. Profiles of the radial velocity dispersions for the DSC compo- 
nent (blue curve), cD component (red curve) and satellite galaxies (symbols 
with errorbars), after stacking all the 44 simulated clusters of our sample. 
Errorbars for galaxies are the 1<t Poissonian uncertainties associated to the 
number of galaxies found within each radial bin. 
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Figure 6. Relation between the DSC stellar mass and cD mass. The solid 
line marks a one to one relation. On average, the DSC has ca. 15% more 
mass than the cD, as indicated by the dashed line. The dotted line is the 
best-fit linear relation, after excluding the outlier corresponding to the cD 
with mass ~ 3 X 1O 12 M . 
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Figure 5. Same as FigurefJ] but comparing the velocity dispersions of DSC 
component, galaxies and dark matter. We also report the value of the nor- 
malization A for the best fit to equation [2] The inlay shows the individual 
velocity histograms for the three components obtained from the most mas- 
sive cluster. 

with a normalization A of 1.58, 2.55 and 17.4 for galaxies, DSC 
stars and cD stars respectively. The value of A we obtain for galax- 
ies is in line with the results presented in lBiviano et 

3.2 The different dynamical tracers 

Although the differences in the velocity dispersion between the 
galaxies and the DSC component as seen in Figure [3] are small - 
basically of the size of the RMS scatter of the correlation - they 
seem to be systematic. Therefore, the question rises of which of 
the two distribution traces the underlying dark matter best. Figure 
[5] shows the velocity dispersion of dark matter, galaxies and the 
DSC component for the full set of analyzed galaxy clusters. Inter- 
estingly, the velocity dispersion of the dark matter falls between 



the one of the DSC component and the one of the galaxies. The 
best fit to equation[2]gives a normalization A of 1.58, 2.06 and 2.55 
for galaxies, dark matter and DSC respectively. The inlay shows the 
velocity histogram (arbitrarily normalized) of the DSC, the galaxies 
and the dark matter within the most massive cluster in the simula- 
tion, clearly revealing the systematic differences among the three 
different components. Note that we had to apply our identification 
scheme, described in section [4] to obtain the velocity distribution 
of the individual stars associated with the DSC. 

3.3 The Mass in the different stellar distributions 

Having the velocity distribution given by our double Maxwellian 
fit, we can also evaluate the total stellar mass associated to the two 
individual populations by just integrating separately over the two 
distribution functions which sums up in the double Maxwellian fit. 
Figure [6] shows the stellar mass associated to the DSC compared 
to the cD component. From this estimate follows that the simula- 
tions predict the mass within the DSC to be, on average, slightly 
larger than the mass in the cD galaxy, by « 15%, as indicated by 
the dashed line in figure [6] From this result follows that the mass 
fraction of the two dynamically identified components associated 
with the DSC and the cD will not depend on the mass of the galaxy 
cluster, being their ratio constant. 

This is different fro m previous findings (e.g. M04, 
Som mer-Larsen et al. (2005)), and it is a feature of the dynami- 
cally based separation of these two components, as we will discuss 
in detail later. 



4 SPLITTING CD AND DSC COMPONENTS 

After detecting two distinct components, which we identify as cD 
and DSC, in the velocity distribution of the stellar population, we 
need a way to assign the individual star particles to one of them. 

To perform this task, we start from a given fiducial radius (ini- 
tially assumed to be a fraction of the virial radius) and, as a first 
step, we calculate the gravitational potential given by all particles 
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Figure 7. As in figure [5] velocity histogram (black line) of the main halo 
stellar component and a double Maxwellian fit to it (gray line) is shown. Red 
and blue histograms show the velocity distribution of the cD and DSC stars 
after the un-binding procedure, together with the two individual Maxwellian 
distributions from the global fit, respectively. 



Figure 9. Relation between DSC and cD stellar mass, after the unbinding. 
The mean relation is the same, but the scatter increases significantly. 
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Figure 8. Main halo (black symbols), cD (red symbols) and DSC (blue 
symbols) stellar density profiles. The red and blue line are a Sersic fit to the 
two distributions, respectively. 



within this fiducial radius and apply the un-binding to the particles 
of the main halo component. As a second step, we separately fit the 
velocity distribution of each of the two stellar populations, classi- 
fied as "bound" and "unbound" in the first step, to a Maxwellian 
distribution, obtaining two velocity dispersions a\, a^ b . 

We compare erf and a\ h to the dispersions obtained from the 
original, double Maxwellian fit performed on the total star compo- 
nent of the main halo, a\ and erf. We then increase or decreases 
the radius for the unbinding procedure, trying to match a\ and erf. 
Using this new radius and re-computing the gravitational potential, 
the procedure is repeated until the velocity distributions of the two 
components converges onto the the two velocity distributions in- 
ferred from the original, global fit: a u t = o\ and Ob — &2- We 
usually adjust the radius to match the value o\ of the DSC compo- 
nent, but if the algorithm does not converge, we try a second time 
using the dispersion 02 of the cD component. The iterative pro- 
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26.0 
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DSC 


516 
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Table 1. Parameters of the Sersic fit for the DSC and cD components 
(columns 1-2) and of the fit to the density profiles for the two components 
obtained from double-Maxwellian fit (see eq. 3). 



cedure stops when a u b/cri (or <Jb/&2) differs from one less than 
a given threshold value. Note that our iterative procedure tries to 
match one velocity dispersion at a time; the fact that, doing so, we 
also obtain the correct value of the other velocity dispersion is not 
given a-priori. Usually both dispersions agree within percent level 
within less than ten iterations. 

Figure [7] shows, similarly to figure [2] the velocity histogram 
of the main halo stellar component (black line) in the same, mas- 
sive cluster of the simulation, fitted by the superposition of two 
Maxwellian distributions (solid lines). We also show the histograms 
of the bound stellar particles belonging to the cD component (red 
histogram) as well as the unbound stellar component associated to 
the DSC (blue histogram). Clearly, the two components are well 
identified and separated by this procedure. 

Once we identify the individual stars in the two distributions, 
we can plot their radial density profile. Figure [8] shows the radial 
density profiles of the cD and DSC component respectively, as well 
as a fitted Sersic profile 



] og (p.)=I -A( 1 f-) 



(V°0 
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We find that both profiles are fitted quite well by individual Sersic 
profiles and therefore the total, radial stellar density profile can be 
described as a superposition of two Sersic profiles. Also, in line 
with previous finding of M04, the profile of the cD component 
is much steeper and less curved than that of the DSC component. 
Therefore, also using our new SUBFIND identification scheme for 
the DSC, we find excess light at large radii with respect to the cD 
Sersic profile. We report in Table[T]the values of the best-fit param- 
eters. 

Comparing the velocity dispersions 0-1,02 inferred from the 
original, double Maxwellian fit to all the main halo stars and the 
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two dispersions at, a u t for cD as well as the DSC we find a 
very tight relation. In almost all cases the iteration was success- 
fully based on the DSC component, which led to a scatter below 
our chosen threshold. Although the dispersions of the inferred cD 
component match very well the corresponding ones of the double 
Maxwellian fit, the scatter is slightly larger. It is difficult to say 
if an improvement of the iteration procedure could further reduce 
this scatter or if it is intrinsically due to our un-binding scheme. In 
effect, the very fact that unbinding star particles using the matter 
contained in a sphere of a given radius divides them into two sets, 
whose velocity distributions match the two component of the over- 
all double Maxwellian distribution is not a-priori guaranteed. We 
regard to this match as a hint that the double Maxwellian shape of 
the overall star particle distribution is really produced by the inter- 
play of two well separated gravitational potentials, the cluster one 
and the central galaxy one. 

Figure|9]shows the relation between the stellar mass within the 
DSC and the cD component, like Figure [6] This time we directly 
use the mass of our two components, as obtained by our unbinding 
scheme. As before, the dashed line indicates an 15% higher mass 
in the DSC than in the cD component. 

A larger scatter in the masses of the cD obtained from the un- 
binding compared to the ones calculated from the velocity distribu- 
tions is clearly visible. Beside this, there is no bias associated with 
our unbinding scheme: it seems able to disentangle two component 
whose properties are identical to those inferred from the velocity 
distribution of the whole stellar population. 

In Fi gure loTTl and l9l we also report results obtained from two 
individual galaxy cluster re-simulations (using different symbols, 
as indicated in the plots), which have 3 and 18 times better mass 
resolution. Such results are perfectly in line with those obtained 
from the cosmological simulations. Thus, our unbinding scheme 
appears to be stable against resolution effects. Note that, at vari- 
ance with other schemes, e.g SKID, we don't have any parameter to 
be tuned when resolution is changed. However, we excluded these 
points from the fits presented in Figure [8] 



5 PROPERTIES OF THE DSC IN THE SIMULATED 
LOCAL UNIVERSE 

We now apply our new un-binding scheme to all 44 massive clus- 
ters in the local universe simulation with at least 20 member galax- 
ies. In all cases, Sersic profiles are a good fit to the radial density 
profiles of cD and DSC populations separately, as already shown 
for one massive cluster in Figure[8] Also, as expected, the cD com- 
ponent in all systems is radially more concentrated than the DSC. 
By inspecting the radial density profiles of the separated compo- 
nents we find that the DSC component in all cases dominates at 
radii larger than 0.03-0.08 R V i r and contributes less than 10% to 
the stellar density in the center. 



5.1 Comparison with previous works 

Fi gure [Tol shows the fraction of stellar mass in the DSC against the 
cluster virial mass. To make a significant comparison with M04 re- 
sults, we also analyzed the local universe simulation using SKID, 
repeating the same procedure described there (see also Appendix 
for details). We remind the reader that M04 analyzed a cosmo- 
logical box of size 192/i _1 Mpc, run with mass and force reso- 
lution similar to those used here, but different physics: in particu- 
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Figure 10. Fraction of stellar mass in DSC as a function of the galaxy cluster 
virial mass. Blue diamonds show the DSC fractions as given by our new 
SUBFIND+un-binding algorithm, red triangles show the DSC fractions as 
obtained using SKID. The red line shows the result of a least-chi-squared 
fit of the SKID fractions to the linear function Fj^gc = b ' Mclus + a > 
which gives a = 0.211, b = 1.0310" 16 , a a = 0.01, a b = 2.5910" 17 , 
X 2 = 0.112. Performing the same fit using the SUBFIND fractions gives 
a value of b compatible with zero at 1.5 a, we thus fit SUBFIND DSC 
fractions to a constant function b = 0.335 ± 0.02 (blue line) 
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Figure 11. The left panel shows the mean formation times of cD (red sym- 
bols) and DSC (blue symbols) as function of the virial mass of the cluster. 
Results from the higher resolution, 3.x and 18x, simulations, are shown with 
triangles and stars, respectively. The right panel shows the histograms of 
the distribution of the mean formation times for the cD and the DSC com- 
ponents, respectively. 

lar, they didn't follow chemical enrichment and didn't use a metal- 
dependent cooling function. 

Nevertheless, we recover the trend of having a larger DSC 
fraction in larger clusters when we define it using SKID, even if 
such a trend is weaker here than in M04. Our new dynamical iden- 
tification of the DSC, however, gives a different results. Although 
the DSC fractions for the most massive cluster agree well, on av- 
erage they are higher than those found using SKID, especially at 
low cluster masses, and, moreover, we can't detect any trend with 
the cluster mass. In Appendix, we show a detailed comparison of 
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Figure 12. Formation time of main sub-halo stars (black histogram), cD 
stars (red histogram) and DSC stars(blue histogram) in the most massive 
cluster in the cosmological simulation. 

one high mass and one low mass cluster analyzed with both meth- 
ods to demonstrate that this effect is indeed related to the fact, that 
for the SKID analysis a scale has to be given a-priori, which could 
lead to this effect. However, we want to stress once more that the 
DSC mass fraction found here should not be regarded as immedi- 
ately comparable to observations. Especially in low-mass clusters, 
the DSC abundance is expected to be dominated by central regions. 
Central regions are exactly those where a dynamical identification 
of a DSC, as that presented here, can most differ from a surface 
brightness based detection of an 1CL components. Dynamical anal- 
ysis of the DSC should instead be regarded as attempts to determine 
its physical properties and its origin. 

5.2 Different stellar populations in cD and DSC 

In line with M04, we find that the stellar populations of the two dif- 
ferent components, however we identify them, differ significantly 
not only in their spatial distribution, but also in their history. 

The left panel of Figure [TJJ shows mean formation times of 
the stars in the cD (red symbols) and DSC (blue symbols) for our 
sample of 44 clusters as obtained using our new procedure. In the 
DSC (and therefore its progenitor galaxies), there is no significant 
star formation within the last 4-5 Gyrs, while the cD population is 
on average younger. The difference in the two age distributions is 
also clear from the right panel of Figure[TT] where we show the his- 
tograms of mean formation times for the cD (red) and DSC (blue). 
This agrees with the finding of M07, that DSC is mainly formed 
during merger processes which also form the central galaxy of 
each cluster. Destroyed satellite galaxies, and pre-processed stars, 
already unbound in merging galaxy groups, goes to the DSC of 
clusters, whereas the cD still can form stars at low rate at recent 
times. In numerical simulation of galaxy clusters, one known issue 
is the excess star formation at the center o f cooling flows, p roduc- 
ing an overly blue central galaxy, see e.g lSaro et al. (2006). Real 
BCGs have a much lower low-z star formation rate, therefore we 
expect this age segregation to be smaller in reality; by how much, 
is stil to be completely understod. Note, however, that this mech- 
anism applies to all mergers at all times. Therefore, also at higher 
redshifts, the merger remnants at the center of proto-clusters will 
have ongoing star formation while, obviously, the mean age of the 
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Figure 13. Fraction of stellar in the main sub-halo, the DSC and the cD 
component, respectively, to the total mass in stars as a function of the frac- 
tion of the stellar mass of the second brightest galaxy in each cluster to the 
cD mass. The results from the high resolution simulations 3.x and 18.x are 
shown with triangles and stars, respectively. 

produced DSC will grow older. This can be seen from Figure Q2] 
where we show, for the most massive cluster in the cosmological 
simulation, the distribution of formation times of star particles for 
all stars in the main sub-halo, cD stars and DSC stars. Star particles 
older than 7 Gyr (i.e., tform < 6) are abundant in the DSC; these 
stars formed at redshifts z > 1 and must have become unbound 
before. At late-times, lack of feedback in the simulations leads to 
ongoing star formation at the center of the cluster gives an excess of 
young stars (tform > 8 Gyr) in the cD. However, even ignoring all 
stars with tf orm > 8 Gyr, the difference between the average for- 
mation time of cD and DSC still remains half of the effect shown 
in figure ITTI 

As before, we added in figure [TT] the results obtained using 
our two re-simulations of individual galaxy clusters. Again, they 
reproduce the trends found in our cosmological simulations, and 
for the same reasons as before, we did not include them in the his- 
tograms. Note that, especially in the very high resolution run, star 
formation is slightly shifted to earlier times, as expected since such 
a simulation resolves the first star-forming structures better. 

We also evaluated mean formation times of cD and DSC stellar 
populations given by SKID for our 44 clusters. The result is similar 
to that shown for SUBFIND in Figure[l2] but mean formation times 
of the two components are closer with the average formation time 
of cDs being 4.9 Gyr and that of DSC, 4. 1 Gyr. Again, this suggests 
that SKID is less efficient in disentangling the two components, and 
assign to some cD star particles which have properties typical of 
the DSC. 

5.3 Origin of the scatter of the DSC 

Both observations and previous numerical works indicated that the 
fraction of the DSC component in clusters having similar mass 
shows a large scatter. Reminding the reader that our cluster set is 
a complete, e.g. volume limited sample, we can try to investigate 
where the scatter in relations of cD and DSC components with the 
cluster mass comes from. Thus, we calculated the the total mass 
for our different stellar components, namely the total amount of 
formed stars, the amount of stars in the main sub-halo and the stel- 
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Figure 14. Same as figure [T"3l but as a function of the formation redshift, 
defined as the redshift at which the mass of the main progenitor is half of 
the mass of the descendant cluster at redshift zero. The results from the 
high resolution simulations 3x and 18x are shown with triangles and stars, 
respectively. 



lar mass of DS C and cD componen ts respectively. In line with pre- 
vious findings l IBorgani et al 1 l2006h . there is a week trend to form 
more stars in smaller clusters, but the scatter of formed stars in 
similar mass systems is quite small. Although the global trends is 
reflected more or less homogeneously in all different components, 
the relative scatter, especially in the main and DSC component is 
quite substantial. 

As seen before, the un-binding procedure induce additional 
scatter by itself; however this does not account for the whole vari- 
ation of the DSC fraction among different systems of similar mass. 
Figure [T3] shows the fraction of all stars in the main sub-halos and 
in the DSC, as a function of the relative importance of the second 
most massive galaxies within each cluster. There is a clear indica- 
tion that the presence of a quite massive second galaxy (indicating 
a young system, where the remain of the last major merger has not 
yet merged with the cD) is related to a less prominent DSC compo- 
nent. Note that this is fully compatible with the low fraction of ICL 
observed in the Virgo cluster (< 10%), which is characterized by 
the presence of two comparable brightest cluster galaxies. There- 
fore the observed Virgo cluster lies even to the right of the shown 
range in figure [T3l 

By constructing the merger tree of our set of clusters we can 
trace back the main progenitor of each of the cluster. Adopting the 
common definition for the redshift of formation of a cluster Zf orm 
(redshift at which the main progenitor of the cluster has half of its 
final mass) we can check how the stellar fractions in the different 
components depend on the age of the cluster. Figure [14] shows the 
stellar fraction of the main sub-halos and the DSC component as a 
function of the formation redshift 2f orm . In old systems, expected 
to be very relaxed, the growth of the cD galaxy seems to be the 
dominant source of increasing the stellar fraction, whereas the DSC 
seems not be evolving significantly. Only in very young systems 
the DSC shows a deficit, most likely because the merging is still 
ongoing and the DSC is not yet released from the progenitors of 
the forming systems. 

Summarizing, when we apply our new un-binding scheme, a 
correlation of the DSC with the cluster mass is not present at all. 
The amount of diffuse stars appears to depend more on the forma- 



tion redshift of the cluster and on its dynamical state: more relaxed 
and evolved clusters have more DSC because they had more time 
to complete the merging of galaxies onto the cD. 



6 CONCLUSIONS 

We showed that the stellar populations of the main sub-halos of 
massive galaxy clusters in cosmological, hydrodynamical simu- 
lations are composed of two, dynamically clearly distinct com- 
ponents. Their velocity distributions can be fitted by a double 
Maxwellian distribution. Including this information into an un- 
binding procedure, we where able to spatially separate the two 
components in a central galaxy (cD) and a diffuse stellar compo- 
nent (DSC). The latter can be associated with the observationally 
very interesting intra cluster light (ICL) component. However, we 
want to stress that there are significant differences in how the ICL 
is measured in observations and how the DSC is inferred in simula- 
tions. Therefore a detailed comparison of the amount of this stellar 
component of galaxy clusters in simulations and observations can 
only be done by mimicking observational strategies, e.g. inferring 
the component from synthetic surface brightness maps. 

On the other hand our un-binding algorithm is ideally suited 
to study the physical processes leading to liberate the DSC as it 
allows to distinguish the DSC and the cD with a minimum set of 
assumptions, namely 

• The stars of the cD component are those gravitationally bound 
to the inner most part of the halo, i.e. to the galaxy itself. 

• The stars of the cD and the DSC reflect two dynamical distinct 
populations. 

These two assumptions naturally reflect our expectations that the 
cD galaxy is a relaxed object sitting in the center of the cluster po- 
tential, while DSC is mainly formed by violent mergers of satellite 
galaxies with the cD, is in equilibrium with the overall gravitational 
potential of the cluster, and still holds memory of the dynamics of 
the satellite galaxies. 

We show that our un-binding scheme is able to disentangle 
two stellar populations, cD and DSC, whose velocity distributions 
separately matches the two components of the double Maxwellian 
velocity distribution which characterizes the stars belonging to the 
main sub-halos of clusters. 

Applying our algorithm to a cosmological, hydrodynamical 
simulation of the local universe including cooling and star forma- 
tion, we found that 

• The velocity dispersion associated to the DSC component is 
comparable to the velocity dispersion of the member galaxies, 
whereas the velocity dispersion associated to the cD is significant 
smaller (by a factor of ~ 3). his implies that the DSC dynamics 
is determined by the general gravitational potential of the cluster, 
while the BCG one feels the local galactic potential. 

• In line with previous findings, the spatial distribution of the 
star particles belonging to the cD are much more concentrated than 
the one of the DSC, which starts to dominate the stellar density at 
radii larger than 0.03-0.08 R v i r - 

• The formation time of the star particles associated with the 
cD is significant smaller than the one of the DSC, indicating that 
- in line with previous findings - the DSC is older (on average 
~ 1.5Gyears). 

• The fraction of the stars within the DSC is only weekly depen- 
dent on the mass of the galaxy cluster, however, we find a larger 
fraction in clusters with a early formation time. This is in line with 
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previous findings that the DSC originates from the late stages of 
merger events of galaxies with the cD. Therefore clusters which 
have had enough time to liberate stars since the last major infall 
show a larger DSC. Interestingly, the mass ratio between the sec- 
ond brightest cluster galaxy and the cD galaxy can be used as a 
proxy for the time till the last major merger. Using such a ratio, we 
obtain an even stronger trend: clusters with a small ratio show a 
large DSC fraction, whereas clusters with a large ratio did not have 
yet enough time to liberate the stars contribution to the DSC. We 
expect that such a trend could also be detected in ICL observations. 

We conclude that the separation of the cD and the DSC in 
simulations, based on our dynamical criteria, is more opportune 
than other current methods, as it depends on less numerous and 
physically more motivated assumptions, and leads to a stable algo- 
rithm helping to understand the different stellar components within 
galaxy clusters. 
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Figure A2. As in Figure[8] the stellar density profile of the cD (red symbols) 
and the DSC (blue symbols) are shown. The red and blue lines represent 
a Sersic fit to the distributions. Diamonds with solid lines and stars with 
dashed lines are the results obtained with SKID and SUBFIND respectively. 
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APPENDIX A: SKID - SUBFIND COMPARISON 

The id entification o f galaxies in M04 and M07 was based on a 
SKID. dStadefeOOll) . For obtaining stable results, SKID parameters 
must be tuned carefully for correctly identifying individual galax- 
ies. Namely, M07 showed that three different smoothing length 
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Figure A3. As in Figure [7] velocity histogram (black histogram) of the 
main halo stellar component and a double Maxwellian fit to it (thin line) 
is shown. Additionally, the red and blue histograms show the velocity dis- 
tribution of the cD and DSC components and the thin lines indicate the 
individual, Maxwellian distributions obtained from the global fit. Light and 
dark colored histograms are the results obtained with SKID and SUBFIND 
respectively. 

must be used, and the results of the three analysis combined, to 
avoid spurious loss of low density galaxies. Moreover, in SKID 
analysis, a typical scale length r must be supplied for performing 
the un-binding procedure. Here (as in M04 and M07), we linked 
such a scale to the gravitational softening e of the simulation: 
r = 3e. M07 showed that this choice is well suited to obtain a 
reasonable separation between the cD and the DSC. 

The differences found between the SUBFIND and SKID for the 
cD vs. DSC separation as a function of mass, could be due to the 
fact that un-binding in clusters of different masses uses the same 
r. In smaller clusters, more DSC particles could be improperly as- 
signed to central galaxies. Such a problem is not present in our new 
SUBFIND scheme. 

One additional small but systematic difference is in previous 
studies based on SKID where always based on all star particles 
within the virial radius, whereas SUBFIND operates on all parti- 
cles belonging to a FoF group. The linking length is usually chosen 
to correspond to the virial over-density in a given cosmology and 
therefore the total mass or the amount of stars is quite similar in 
both analysis, but the FoF group usually is elliptically elongated 
and therefore the particles sets the two algorithms start from are 
not exactly the same in the outer region of the cluster. 

Al A detailed case study for a massive cluster 

Here we compare the analysis of the high resolution simulation 
(e.g. the 3x run) applying both methods to check the differences 
details when applying the two algorithm. 

Figure lATl shows a direct comparison of all star particles in the 
different components as identified by SKID (left column) and SUB- 
FIND (right column) within half of the virial radius of the cluster. It 
is quite encouraging that practically all galaxies have counterparts 
in both methods. In general, it seems that the un-binding applied 
by SUBFIND tend to unbind a little bit more the outer envelope of 
satellite galaxies, which in return appear as small over-densities in 
the DSC component. On the contrary, the cD identified by SUB- 
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Figure Al. All star particles within R v i r /2 of the high resolution (Jjr)cluster belonging to individual galaxies (black dots in the upper row), the cD component 
(red dots in the upper row) and the DSC component (lower row). The left column is obtained applying SKID, the right column is our SUBFIND result. 



FIND has smoother boundaries, while it looks more truncated in 
the SKID analysis. Nevertheless, the mass of the cD component is 
very similar, namely 1.28 x 1O 13 M /h versus 1.21 x 1O 13 M /h 
in SKID and SUBFIND respectively. Note that the cD velocity dis- 
persion is similar in both analysis. SKID truncates the galaxy sim- 
ply because its farther stars are not assigned to the cD, but to some 
satellite galaxy to which they don't belong, during SKID FoF phase. 
Then, they become unbound in the un-binding phase. SUBFIND as- 
sign such stars to the main sub-halo, and lately, they are not unbind 
from the cD. 

The DSC looks extremely similar in both algorithms, which is 
also reflected in a quasi indistinguishable density profile as shown 



in figure IA2I Here also the cD component shows excellent agree- 
ment in the inner part, again the differences in the outer parts are 
visible. We note that also in the outer parts the DSC components 
in both methods fall on top of each other, despite the previously 
discussed small differences in the underlying particle sets. 

Figure lA3l shows the velocity distribution in the two compo- 
nents. Already in the SKID analysis the difference in the dynamics 
of the two components are visible. Although again very similar, the 
SUBFIND results are more closely following the decomposition of 
the main stellar component into two Maxwellian distributions, as 
expected. Again, despite the different underlying algorithms, the 
total velocity distributions of all star particles used in SKID and 
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Figure A5. As in Figure |A3| we show the velocity histogram (black) of 
the main halo stellar component and a double Maxwellian fit to it (thin 
line). Additionally, the red and blue histograms show the velocity distribu- 
tion of the cD and DSC components and the thin lines indicate the individ- 
ual, Maxwellian distributions obtained from the global fit. Light and dark 
colored histograms are the results obtained with SKID and SUBFIND re- 
spectively. The thin lines are individual Maxwellian fits to the cD and DSC 
components identified by SKID. 

SUBFIND is quasi indistinguishable and also leads to the same dou- 
ble Maxwellian fit. Therefore, we plotted in figure lA3l onlv the total 
distribution and the Maxwellian fits to SUBFIND data to avoid con- 
fusions. 

A2 An example of a low mass cluster 

Analyzing a smaller cluster (M„, r ~ 2 x 10 14 M@ / h) from our 
simulations reveals the differences between the two algorithms for 
such less massive systems. Although the individual satellite galax- 
ies match as nicely as before, SKID assigns much more stars to 
the central cD than SUBFIND leading to the contrary behavior for 
the DSC component, which is clearly visible in figure IA4I which, 
similarly to figure [ATI shows a particle by particle comparison be- 
tween the two algorithms. The differences in disentangling the cD 
and DSC component is also clear from figure IA5I which shows 
the velocity distribution histograms of the inferred cD and DSC 
components. Whereas the SUBFIND results follow nicely, by con- 
struction, the individual components of the double Maxwellian fit 
to the global velocity distribution, the SKID results deviate signif- 
icantly. In addition, we fitted individual Maxwellian distributions 
to the SKID results, shown as thin lines. It appears quite obvious 
that even such fits are not as good as the ones obtained with our 
new SUBFIND algorithm. Especially the cD component identified 
by SKID shows a significant tail towards higher velocities which 
can't be represented by a single Maxwellian distribution, indicat- 
ing that the star particles associated to the cD by SKID are not a 
single, dynamical component. 

To obtain better result with SKID, one should probably tune 
the r parameter on a cluster-by-cluster basis. Such tuning is not 
needed with our new SUBFIND scheme. 
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Figure A4. Same as figure lATl but for the small cluster where we found significant differences between the cD and DSC separation when applying SKID (left 
column) and SUBFIND (right column) respectively. 



